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Abstract 

A generalized macroscopic nonlocal theory of sound propagation in rigid-framed porous media saturated with a 
viscothermal fluid has been recently proposed, which takes into account both temporal and spatial dispersion. 
Here, we consider applying this theory capable to describe resonance effects, to the case of sound propagation 
through an array of Helmholtz resonators whose unusual metamaterial properties such as negative bulk moduli!, 
have been experimentally demonstrated. Three different calculations are performed, validating the results of the 
nonlocal theory, relating to the frequency-dependent Bloch wavenumber and bulk modulus of the first normal 
mode, for ID propagation in 2D or 3D periodic structures. 

Resume 

Description nonlocale de la propagation du son dans une chaine de resonateurs de Helmholtz. Une 

theorie macroscopique nonlocale generate de la propagation du son dans les milieux poreux a structure rigide 
satures par un fluide viscothermique a recemment vu le jour. Tenant un compte complet des dispersions tant 
temporelles que spatiales, elle decrit entierement les resonances. Ici, nous I’appliquons au cas de la propagation 
du son dans un reseau de resonateurs de Helmholtz, dont les proprietes non usuelles (modules de compressibilite 
negatifs) ont ete etablies experimentalement. Trois calculs differents sont presentes, qui valident les resultats de 
la theorie nonlocale, relatifs au nombre d’onde et module de compressibilite fonctions de la frequence, du mode 
de Bloch principal (le moins attenue), pour une propagation ID en geometries periodiques 2D ou 3D. 


Key words: Helmholtz resonators ; Acoustic metamaterials ; Nonlocal description; Spatial dispersion; Viscothermal fluid; 
Negative modulus 

Mots-cles : Rsonateur d’Helmholtz; Metamatriaux acoustiques ; Description nonlocale; Dispersion spatiale; Fluide 
viscothermique; Module de compressibilit ngatif 


Email addresses: nnematiSmit.edu (Navid Nemati), akumrSmit.edu (Anshuman Kumar), 
deiiis.lafargeOuniv-lemans.fr (Denis Lafarge), nicfangOmit.edu (Nicholas X. Fang). 


Preprint submitted to Elsevier Science 


March 18, 2015 





1. Introduction 


We employ here a generalized macroscopic nonlocal theory of sound propagation in rigid-framed porous 
media saturated with a viscothermal fluid [1] to describe the behavior of an acoustic metamaterial made 
of an array of Helmholtz resonators filled with air (see Fig. 1 left). Inspired by the electromagnetic theory 
and a thermodynamic consideration relating to the concept of acoustic part of energy current density, this 
theory allows to go beyond the limits of the classical local theory and within the limits of linear theory, 
to take into account not only temporal dispersion, but also spatial dispersion. In the framework of the 
new approach, an homogenization procedure is proposed, through solving two independent microscopic 
action-response problems each of which related to the effective density and effective bulk modulus of the 
material. Contrary to the classical (two-scale asymptotic) method of homogenization, there is no length- 
constraint to be considered alongside of the development of the new method, thus, there would be no 
frequency limit for the medium effective properties to be valid, and also materials with different length 
scales can be treated. The homogenization procedure offers a systematic way of obtaining the effective 
properties of the materials, regardless of their geometries. These characters of the nonlocal approach give 
the possibility to describe the porous media with specific geometries causing metamaterial behavior. A 
metamaterial with periodic structure will be studied: two-dimensional and three dimensional chain of 
Helmholtz resonators connected in series. 

By the local theory we refer to space locality. Nonlocality in time, or temporal dispersion, has been 
already taken into account through models for wave propagation in porous media [2,3,4,5]. That means, 
in Fourier space the effective density and bulk modulus depend on the frequency w. In other terms, the 
field dynamics at one location retains a memory of the field values at this location but is not affected by 
the neighboring values. The local description is usually based on retaining only the leading order terms in 
the two-scale homogenization method [6,7,8,9,10,11,5] using an asymptotic two-scale approach in terms 
of a characteristic length of the medium, the period L in periodic media, which is supposed to be much 
smaller than the wavelength A [12,13]. Efforts have been performed to extend the asymptotic method of 
homogenization to higher frequencies for the periodic composite materials [14,15] and rigid porous media 
[16] by introducing another type of scale separation to which the asymptotic multi-scale procedure applies. 
Enhanced asymptotic method has been adapted to describe sound propagation in rigid porous media with 
embedded damped Helmholtz resonators [17] exhibiting scattering different from Bragg scattering at high 
frequency in periodic media. 

An effective medium approach has been proposed for periodic elastic composites based on surface 
responses of a structural unit of the material [18], which can describe the macroscopic parameters beyond 
the frequencies within the long wavelength limit. Unlike the classical methods, based on the introduction of 
two-scale asymptotic expansions, or coherent potential approximation [19] based on the effective-medium 
parameters minimizing scatterings in the long-wavelength limit, the homogenization scheme presented 
in [18] uses matching the lowest-order scattering amplitudes arising from a periodic unit cell of the 
metamaterial with that of a homogenized material. As such, local resonant scattering can be captured 
as well by the latter in the elastic metamaterials. Enhanced asymptotic method of homogenization has 
been developed to provide the weak nonlocal effects as a small correction to the local behavior [20]. 
An approach has been presented [21] for random elastic composites based on ensemble averaging of 
the material responses to a body force giving rise to effective parameters of the medium depending on 
frequency and wavenumber. By this method the case of periodic media can be treated as well. 

The nonlocal theory we use here, takes fully the temporal dispersion and spatial dispersion into ac¬ 
count. The medium is assumed unbounded and homogeneous in the stationnary random statistical sense, 
therefore, the spatial dispersion refers only to the dependence of the permittivities, i.e. effective density 
and bulk modulus, on the Fourier wavenumbers k present in the macroscopic fields [22]. The theory 
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applies with some care to a periodic medium; in particular it gives the Bloch wavenumbers and defines 
Bloch impedances: indeed, it applies exactly to the ensemble defined by the random translation of one 
periodic sample. The materials susceptible to show the nonlocal behavior may be classified into two main 
groups regarding their microgeometry. The first comprises the materials which exhibit this behavior in 
sufficiently high frequency regime. The second one concerns materials with microgeometry constituting 
the resonators, which exhibit spatial dispersion phenomena even at not very high frequencies; the reso¬ 
nance phenomena act as a source generating nonlocal behavior. In this article we investigate the second 
type of these geometries in the form of daisy chained Helmholtz resonators. We will see the first one in a 
forthcoming paper, where ID propagation in a two-dimensional lattice of rigid cylinders will be studied. A 
material made of an array of Helmholtz resonators filled by water has been studied experimentally, which 
has been found to show negative bulk modulus in the resonance frequency range [23]. Later, Helmhotz 
resonators as structural units were used to designe novel metamaterials for focusing ultrasound waves [24] 
and broadband acoustic cloaking [25]. 

Here, we will apply the nonlocal theory to quantatively describe the macroscopic dynamics of such 
a metamaterial filled with air as a viscothermal fluid, in 2D as well as in 3D. For the 2D case, using 
a simplified analytical solution of the complete equations, we will present how to obtain the nonlocal 
efffective density and effective bulk modulus. When these effective parameters satisfy the dispersion 
equation based on the nonlocal theory, we can compute the wavenumber of the least attenuated mode, 
among other modes. We can then check that the wavenumber resulting from the macroscopic nonlocal 
theory coincides with the Bloch wavenumber propagating and attenuating in the medium. The Bloch 
solution is determined using the same simplifying solving as in the nonlocal modeling, thus the results 
based on the two calculations should be comparable. Finally, as a check of the validity of the simplifying 
assumptions introduced in our modeling calculations, we have performed direct Finite Element Method 
(FEM) computations based on the exact equations in the framework of nonlocal homogenization. 

In section 2, we review briefly the general framework of the nonlocal theory which will be used after¬ 
wards. The microscopic equations governing sound propagation in a rigid porous medium will be recalled, 
before mentioning the macroscopic Maxwellian equations describing the macroscopic nonlocal dynamics 
of the homogenized equivalent fluid. In section 3, we will see the nonlocal modeling allowing the calcu¬ 
lation of the effective parameters and the wavenumber of the least attenuated wave in the medium. The 
direct calculation of the Bloch wavenumber, using the similar simplifications, is presented in section 4. 
Section 5 is devoted to the results of the three different calculations in 2D, and also the results based on 
the nonlocal modelling and Bloch wave calculations which have been generalized to 3D structures. 


2. Genereral framework of the nonlocal theory 

In the following, we state the microscopic equations applied at the pore level, and the nonlocal 
Maxwellian macroscopic equations that describe the dynamics of the material as a homogeneous equiva¬ 
lent fluid medium. Then, we recall briefly the upscaling procedures allowing to obtain the frequency and 
wavenumber dependent effective parameters of the macroscopic equivalent fluid medium, i.e. effective 
density and effective bulk modulus [1]. 

2.1. Microscopic equations 

The dynamics of a small amplitude perturbation in a rigid-framed porous material filled with a vis¬ 
cothermal fluid is governed by the linearized equations of the mass, momentum, and energy balance, and 
a general fluid state equation as follows: in the fluid region 


3 




"I “ 



Figure 1. Left: illustration of a 2D array of Helmholtz resonators. Right: a periodic cell of the structure, with L = 1cm, 
E = 0.2L, cr = 0.015L, and / = 0.15L. 
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with boundary conditions 


= 0, T = 0 (2) 

applied to the fluid/solid interface 9V, where v, b = p/po, p and t, are the fluid velocity, excess con¬ 
densation, thermodynamic excess pressure, excess temperature, respectively, and p is the excess density. 
The fluid constants po, 77 , C, 7, xo^ Po, Cp^ To, n, represent the ambient density, first viscosity, second vis¬ 
cosity, ratio of heat coefhcients at constant pressure to constant volume Cpjcy, adiabatic compressibility, 
coefhcient of thermal expansion, specific heat coefficient at constant pressure, ambient temperature, and 
coefficient of thermal conduction, respectively. 


2.2. Macroscopic Maxwellian acoustics 


Before going through the macroscopic equations for sound propagation in rigid-framed porous media, 
and the homogenization procedure, we will precise the notion of field averaging in the nonlocal approach. 

Averaging: the present macroscopic theory is statistical in nature and has been developed in principle 
for fluid-saturated rigid-framed media which are homogeneous in an ensemble-averaged sense; this is the 
case of stationnary random media. The macroscopic properties represented in the theory refer to the 
ensemble of realizations. Thus for example, the propagation constants of the medium would refer to the 
propagation constant of coherent waves in multiple-scattering theory. Here, the material we wish to study 
is not defined by stationary random realizations. It belongs on the contrary to the important class of 
periodic materials. The macroscopic theory can still be applied, however, if we consider the ensemble 
obtained by random translation of one sample. It turns out that the ensemble-average () properties of 
the space are, in this case, precisely computable by spatial averaging over a periodic cell in a single 
realization. This, in a sense, reminds of ergodicity in the stationary random case. 
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The macroscopic condensation and velocity are defined as the average of pore scale microscopic fields: 
V = (v), and B = (b); average over the periodic cell in the case of the periodic media. The electromagnetic 
analogy suggests that the system of macroscopic equations can be carried through by introducing new 
Maxwellian fields H and D, and also operators p and such that 


Field equations: 


dB 


v-y = 0, 


dD 


= -VH 


( 3 ) 


Constitutive relations: D = pV, H = X~^B 

where the integral operators of density p and bulk modulus x~^ are defined by 

i-t 


D{t,r) = f dt' f dr'p{t — t',r — r')V{t',r') 

J —OO J 

H{t,r) = f dt' f dr'x~^{t — t',r — r')B{t',r') 

J —OO J 


( 4 ) 


(5a) 

(5b) 


We notice that the kernels p and x~^ depend on the difference t — t' and r — r', which is due to the 
homogeneity in time and material space. That is why we can write (5a) and (5b) in the Fourier space, 
respectively, as 

D{LO,k) = p{uj,k)V{uj,k), H{LO,k) = x~^{i^,k)B{uj,k) (6) 

In nonlocal theory, the macroscopic H field is defined through the Poynting-Schoch condition of acoustic 
part of energy eurrent density [1,26]: 

S = HV = (pv) (7) 

Regarding Eqs.(5) and (6), it is visible that the theory allows for both temporal dispersion, shown by 
integration over time variable t' in physical space and frequency dependence in Fourier space, and spatial 
dispersion, shown by integration over space coordinates r' and wavenumber dependence in Fourier space. 
We will recognize the quantities in physical space (t, r) and Fourier space (w, k) by their arguments. Now, 
in order to clarify the relationship between constitutive operators and microgeometry, the kernel functions 
p{u),k) and x“^(a;,fc) are needed to be determined, by introducing the procedures coarse-graining the 
dissipative fluid dynamics of the pore scale. 


2.3. Procedures to compute effective density and bulk modulus 

In the ID case of macroscopic propagation along a symmetry axis, for instance x-axis with the unit 
vector i, we will have D = Dx and V = Vx^ r = xx, and k = kx in the above equations (3-7). To 
determine the Fourier functions p{uj, k) and x“^(w, k) for the ID acoustic propagation in a medium with 
porosity (j), we solve two independent action-response problems. For computing the effective density we 
consider the macroscopic response of the fluid subject to a single-component (w, k) Fourier bulk force. 
The effective bulk modulus is related to the response of the fluid subject to a single-component Fourier 
rate of heat supply. 

Two sets of equations to be solved: the two systems of equations to be solved are written as 

In the fluid region V/: 


5 



(8a) 
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dr o rr. dp 2 , 

Pocp^ = PoTo-^ + kV T + 


Fe 




Added for determination of density 


Added for determination of bulk modulus 


IXoP = b + /3 ot 


(8b) 

(8c) 

(8d) 


On the fluid/solid interface dV: 

V = 0, T = 0 (9) 

For convenience the excitation amplitudes are written as: Qe-^‘^t+ikx _ f]gTo{d/dt) and 

p^-iut+ikx _ ^p^-zuit+tkxy Here, it is important to note that the excitation variables oj and k are 
set as independent variables. The solutions to the above systems for the fields p, b, r, and components of v 
take the form p{t, r) = p{uj, k, and so on. Recall that the theory is formulated for a geometry 

that is stationary random, and the averaging ( ) is the ensemble average. Thus here, the amplitude fields 
v{uj, k, r), p{uj, k, r), b{u}, k, r), and t{uj, k, r), are stationary random functions of r. Passing to the case of 
periodic geometry, we can limit ourselves to considering one periodic sample. The fields become periodic 
functions over a cell, and ( ) is interpreted as a volume average over a cell. 

Effective density and bulk modulus: once the two systems of equations are solved independently, 
using the right hand Maxwellian macroscopic equations in (3) and (4), we arrive at the following expres¬ 
sions for the nonlocal effective density and bulk modulus 


p{uj,k) = 


k{r + p{uj,k)) 
uj{v{uj,k,r)) 

^ PiLO,k)+V 
{b{u},k,r)) + (j)-fXoP 


(10a) 

(10b) 


where P{v) = (pv), which has been inspired by (7). 

Wavenumbers: contrary to the case of local theory, here, since we fully take into account spatial 
dispersion, several normal mode solutions might exist, with fields varying as Each solution 

should satisfy the following dispersion equation 

p(w,fc)x(w,fc)u;^ = (11) 

which is easily extracted from the Maxwellian macroscopic equations. With each frequency w, several 
complex wavenumbers ki{ui), ^(ki) > 0, I = 1 , 2 ,..., may be associated. 

In what follows, with the aim of obtaining the nonlocal effective density, effective bulk modulus, and 
wavenumber of the least attenuated mode, we will apply this theoretical framework in analytical simplified 
manner, to a 2D array of Helmholtz resonators, illustrated in Fig. 1 right, exhibiting resonance phenomena 
which result in metamaterial behavior. 


3. Nonlocal modeling for 2D structure 

We proceed to determine the functions p(uj, k) and x“^(a;, k) sufficiently precise to give an appropriate 
modeling of the least attenuated mode, which results then in purely frequency dependent functions p{uj) 
and x“^(a;). To this aim, we need not consider in full detail the microscopic fields v and p. In the 
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waveguide t and cavity c, instead of the microscopic fields, we can use the mean values Vt(c) = {v)s ■ x 
and Pt{c) = {p)s, where ( )s denotes the average at a given x over the waveguide or the cavity width; 
and in the neck n, we can use the mean values Vn = {v)s ■ y and Pn = {p)s^ where ( )s denotes the 
average at a given y over the neck width, and y is the unit vector in the y direction. At the same time, 
we make some simplifications consistent with describing the propagation of these averaged quantities in 
terms of the Zwikker and Kosten densities p{uj) and bulk modulii in the different slit portions. 

These depend only on the slit half-widths, which we shall denote by s*, s„, and Sc, in the tube, neck, 
and cavity. The different slit-like tube portions are illustrated in Fig.2. The main tube t is divided in two 
Zwikker and Kosten ducts, a left duct, and a right duct, oriented in the x direction. The same separation 
is made for the cavity c, whereas the neck n is not divided but seen as one Zwikker and Kosten duct 
oriented in y direction. 


3.1. Determination of nonlocal effeetive density 


Considering the periodic cell of Fig.l right, and the corresponding cell average operation ( ), we look 
for the response of the fluid when a harmonic driving force f{t,x) = in the direction of is 

applied. If we can determine the microscopic response velocity and pressure fields u, p, then we will have 
the function p{uj,k) through the relation (see Eq.(lOa)) 


p{uj,k) 


f - ffcP(a;,k) 
—iui(v{uj, k, r)) 


( 12 ) 


with P(w, k) = (pv) / (v), where the v is the a;-component of the microscopic velocity v. 

In [26, Appendix], the Zwikker and Kosten local theory is expressed for tubes of circular cross-section. 
For 2D slits, exactly the same general principles of modeling may be used; only some details of the 
calculations are changed. In particular, the Bessel functions Jq and Ji are replaced by cosh and sinh 
functions. Zwikker and Kosten’s effective densities Pa{oj) and bulk modulii x~^(ijj) in the guide, neck and 
cavity, will be [27] 



Figure 2. Illustration of slit portions and plane waves propagating in different parts of the resonator. Different positions are 
indicated by m, and different amplitudes by Am-, = 1, •••, 10. 
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Pa{uj) = po 


1 - 


tanh iwpoSq 
^-iuposl/p 


Xa = 7^0 


1 + (7- 1) 


tanh —iojpQCpsI^/ 
yJ-iupoCpsl^jn 


-, -1 


(13) 


for a = t,n,c, where the indexes t, n, and c are related to the tube, neck, and cavity respectively, po, 
and Pq the fluid pressure at rest. The corresponding wavenumbers ka(to) and characteristic admittances 
Ya{uj) are expressed as ka = iojca-, and Ya{ui) = 2sq,/(pqCq), for a = t,n,c, where Ca = ^ly/PaXai is 
the corresponding Zwikker and Kosten’s phase velocity. Notice that we include the slit width 2 sq, (resp. 
S, cr, and L — S — 2/ in the resonator, see Fig. 1, left) in the definition of the characteristic admittance, 
because it simplifies the subsequent writing of continuity conditions. 

We start writing the Zwikker and Kosten’s equations in the different parts of the periodic cell. For the 
tube and the cavity, i.e. , a = t, c, we have 


■ Pa(w) dP, 

——-—\4v = — 


5a 


dx 


/e 


ikx 


lUbaXa{^)Pa = 


(14a) 

(14b) 


where, Va = V^Sa is the flow rate field across the cross section 5^, with I 4 the x-component of the 
velocity in the sense of Zwikker and Kosten (averaged over the section), and Pa is the Zwikker and 
Kosten’s pressure. In the neck, the external excitation having no y-component 


. Pn(w)_ 

lUJ - 14 


iuiaXn{u})Pn 


dPn 

dy 

dy 


(15a) 

(15b) 


where, = V^cr is the flow rate, with Vy the y-component of the velocity, and P„ is the Zwikker and 
Kosten’s pressure in the neck. 

The general solution of the non homogeneous equations in the tube and the cavity, {Pa, Va), a = t,c, 
is written as the sum of the general solution {Pa,h, Va,h) of the homogeneous equations and a particular 
solution {Pa,p, Va,p) of the non homogeneous equations. A general solution of the homogeneous equations 
(14) is written as 


/ Pa./^ j ^ I M (16) 

\Va,h j \Ya j \-Ya j 

where A+ and A~ are the amplitudes of the plane waves in direction of the positive x-axis and negative 
x-axis, respectively. The following particular solution can be considered 



(17) 


where Ba and Ca are four constants (for each w) to be determined. Substituting (17) in (14) gives 
the four constants Bt = ik/{oj'^ptXt — k'^), Ct = ioJXtP‘/{^"^PtXt — k"^), Be = ik/{uj'^pcXc — k'^), and 
Cc = iu!Xc{L — S — 21 )/{uApeXe — k'^)- The particular solution is the same in the left and right portions of 
the tube and the cavity. On the contrary and because of the presence of the neck, the general solution will 
have different amplitude constants in the left and right portions. Thus, the general solution of Eqs.(14) 
can be written as 




(18a) 

(18b) 


where (18a) with amplitudes Ai and corresponds to the left part of the tube, and with amplitudes A^ 
and A 4 to the right part (Fig. 2); similarly for (18b): A 7 and Ag for the left part of the cavity, and Ag 
and ^10 for the right part (Fig. 2). These eight amplitudes are to be determined. The general solution of 
Eqs.(15), {Pn,Vn) has the form 



where Ag and Aq are the neck amplitude-relating constants to be determined (Fig. 2). 

Indeed, in the framework of our simple plane-wave modeling, there are 10 relations concerning the flow 
rate and pressure, which are assumed to be verified. These continuity relations involve the values of the 
fields at different locations indicated by numbers m = 1,..., 10, in Fig.2. We now proceed to write them. 

The Bloch condition results in ^ and Then ^ 

-h ^26*'=*'^/^) and We assume 

the continuity of the pressure at the junction (2)-(3), P^ = P^ \ then Ag + A4 = 2I1 -|- Ag. We 
assume the continuity of the pressure at the junction (5)-(2), Pn^'^ = Pf^\ then Age~''^''^l'^ -I- = 

Ai+Ag + Bt- The flow rate at the junction (2)-(3)-(5) is assumed to verify — = Vn^\ which yields 

Yt (Ai — yl2 — ^3 -I- A'^^ = Yn . The continuity of the pressure at the junction (6)- 

(7), Pn'^ = Pc'^'’ results in _|_^gg-*fenV2 — a^j AgY Be- The flow rate at the junction (6)-(7)-( 8 ) 

is assumed to verify Vn^'^ + Vc'^^ = Vc^\ Yn -|- ^(217 — ^s) = lc( 7 l 9 — Aig). The 

pressure is continuous at (7)-(8), P^'^ = then, Af + Ag = Ag Y ^lo- Tbe flow rate vanishes at 
the interface solid-fluid, Vc^ = 0 , we have Yc [7l7e“*^‘=^^“b/2 _ ^gg*fec(i-0/2j = _(j^^-ik{L-i)/2 ^ 

flow rate vanishes at the interface solid fluid, = 0 , we have W [7l9e®^‘=^^“b/2 _ ^^gg-*fcc(i-0/2j _ 

_( 7 ^gife(-£'-b/ 2 ^ 

As such, we have 10 equations for 10 unknown amplitudes Ai,.., Aiq. Once these are determined, we 
will have all the Zwikker and Kosten’s fields through Eqs.(18a), (19), and (18b). At this point, we can 
easily obtain the cell averages (v) and (pv). Let us start with {v) regarding the fact that the Zwikker and 
Kosten’s flow rate has no component along the x-axis 

1 / gO fL/2 M p{L-l)/2 \ 

(v) = — I / Vtdx+ Vtdx+ Vcdx+ Vc dx \ (20) 

^ \J-L/2 Jo J-(L-i)/2 Jo J 

Similarly, we can compute {pv) through the following relation 

1 / gO |•L /2 pO |■{L-l )/2 \ 

(pv) = / PtVtdxY / PtVt dx + / PcVcdxY / PeVedx] (21) 

P \J-L/2 Jo J-(L-l)/2 Jo J 

Now, we can obtain explicitly the effective density function p(w, k) through Eq.(12). In the next section, 
the effective bulk modulus is computed in a similar way but with a different excitation term, and with 
exactly the same conditions on the flow rate and pressure fields at different junctions. 
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3.2. Determination of nonlocal effective bulk modulus 


Considering the periodic cell (Fig.2), when a harmonic heating Q{t, x) = = —icoPfTfP 

is applied in the medium, we write the Zwikker and Kosten’s equations, in each part of the resonator: 
tube, neck, and cavity. The aim is to obtain the function k) as it is indicated in Eq.(lOb). In the 

main tube and the cavity, for a = t,c, we write 


. Pa(w) dP, 

— ICO —- Va = — 




dx 


iu>SaXai>^)Pa + i^jSa (Xa(w) “ 7 Xo) P = 


dVc, 

dx 


( 22 a) 

( 22 b) 


The second term in the second equation might not seem to be obvious but follows the very procedure of 
obtaining (10b). In the neck, the equations are written as 


. p„(w) dP, 
lu) -14, = 


dy 


rwcrx„(w)P„ + iiva (x„(w) - 7x0) P ^ = 


dVn 

dy 


(23a) 

(23b) 


where the term V comes from the averaging of Q over the neck section. Here also, the second 

equation might not appear obvious, but follows the procedure of the determination ( 10 b) seen in nonlocal 
theory. 

As in the section 3.1, the general solution of the non homogeneous equations ( 22 ) in the right or left part 
of the tube and the cavity, is written as the sum of the general solution {Pa,h, Va,h) of the homogeneous 
equations and a particular solution {Pa,p, Va^p) of the non homogeneous equations. A general solution of 
the homogeneous equations (22) is written as Eq.(16). The following particular solution can be considered 



(24) 


where Ba and Ca are constants to be determined. Substituting (24) in (22) gives the four constants Bt = 
- 7Xo)/(fc^ - w^ptxt), Ct = ookixt - 7Xo)S/(fc2 _ = uj'^PciXc - 7Xo )/iP - w^PcXc), 

and Cc = (jjk{xt — 7Xo)(A — E — 2L)/{k'^ — oj'^pcXc)- Thus, the general solution Eqs.(22) can be written 
as Eqs.(18), replacing / with V. The amplitudes Ai, A 2 , A 3 , A 4 , A 7 , Ag, Ag, and Aig (Fig. 2) are to be 
determined. 

As for the tube and the cavity, the general solution of the non homogeneous equations (23) in the neck, 
is written as the sum of the general solution {Pn,h, Vn,h) of the homogeneous equations and a particular 
solution {Pn,p, yn,p) of the non homogeneous equations. We can find a particular solution in the following 
form 



(25) 


where Bn and Cn are two constants which will be determined by substituting (25) in (23): 

Bn = (fljku') {"iXoIXn — 1) sin(fcCT/2), and Cn = 0. To obtain the above expression for H„, the average 
{P'kx'^^ has been easily calculated 



Akx 


dx = -— sin 
ka 


ka 

T 
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Thus, the general solution of Eq.(23) in the neck can be written as 










(26) 


where and Aq are amplitude-relating constants to be determined (Fig. 2). 

As in the previous section 3.1, in the framework of our modeling, there are 10 relations which are 
assumed to be verified, allowing to relate the flow rate and pressures at different indicated points in 
Fig.2. These relations result in 10 equations by which we can compute the amplitudes Ai, ..., Aiq. 
Consequently, all Zwikker and Kosten’s fields will be found. The averages (v) and (pv) are found through 
rewriting the equations (20) and (21) for the actual fields. We need also the expression for (b) to obtain 
We have 


—iu} {b) = —^ ( V • u dxdy 


1 

lA 

IL 

'T2 


v-ndS = -^ ( 




ktL 


k+L / -k+L -ktL .k+L ■ktL\ 

2iCt sin (-Aie"*— + Aae*— + Age*— - ^ 46 -*—) 


where n is the normal unit vector outward from the border of integration. 

Now, we can obtain explicitly the effective bulk modulus function x“^(u;,fc) through Eq.(lOb). 


4. Bloch wave modeling 

In this section we directly seek, without using the principles of the nonlocal macroscopic theory but 
within the same plane wave modeling, the macroscopic Bloch wavenumber ks oi the least attenuated 
wave propagating in the direction of positive x-axis, such that 





(27) 


To the field constituted of 10 Zwikker and Kosten’s slit waves, as illustrated in Fig.2, are associated 10 
complex amplitudes Ai,..., Aiq. As before, on these 10 amplitudes there are 2 relations (27) expressing 
the Bloch condition, and 8 relations expressing the continuity equations. All these relations are now 
homogeneous relations, so that nontrivial solutions will be obtained only if the determinant vanishes. 
This condition will give the Bloch wavenumber ks- 

The first step is to determine the entrance admittance of the resonator W = Vn^'^ / Pn'^. The general 
solution of the homogeneous form of Eqs.(14) for the cavity, a = c, without the forcing term, is written 
as 



where A 7 and Ag are the amplitudes of the waves in the left part of the cavity, and Ag, Aio are the 

f'7\ /Q\ 

amplitudes of the waves in the right part. Regarding the above equation, the three conditions Pc —Pc , 
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= 0, and = 0, result in the three following relations As = Are Ag = Are 

and Aio = Ay. Using (28) there follows pP = A 7 (l - , vP = Y^Ay (l + , and 

— !)■ Then, we can obtain the expressions for Pn^^ and Vn^\ through already 
indicated continuity conditions Pn'^ = Pi’^\ and which, subsequently, yields the 

impedance Yq = Vn^'^ / Pn'^ 

_ p-ikciL-l) 

^6 = ( 29 ) 


1 -)- Q-ikc{L-l) 


Once and Ui®^ are known, we can obtain pl®^ and through 


cos knl -—sinfc„A /p^®) 
—iF„sinfc„,/ cos knl / 


Thus, the impedance of the resonator Yr is expressed as 


-iF„ sin A:„Z + Ue cos 

Ir — Y 

COS knl — i— sin knl 

Now, we look for the macroscopic wavenumber ks- The following relations are satisfied in the right and 
left part of the tube 

(pjw)\ ( cM Pyw-b 

I I “ Y ■ hL ' k,L I J,.(2),(4) I 

\Vt J \-iYtsm— cos— J ykt J 
Making use of Eq.(27), the above equations result in 


_ ^iksL 


cosktL -^sinfctP\ fPj;^'' 
—iYtsinktL cosktL j \A 


On the other hand, as we have seen before, the three following conditions are assumed in the res¬ 
onator: P(®^ = P/^\ Pn'^ = Pt‘^\ and = Vn^\ We have immediately p/®^ = = 

(1/W) {yi^^ — Writing the two equations resulting from (33), and eliminating P^^®^ and P^'^ in 

these equations, gives 


1 _ gjfcBi I ^ _ A sinfctL) (1 - cosktL) 

Yr \Yr Yt J Yr 


^ sin/ctP — cosfctPj 1 — sinfctP j \^t / \^/ 

The determinant of the coefficient matrix must vanish, if the above equations have non-zero solutions. This 
yields a second degree algebraic equation -1-1 = 0, with D = 2 cos ktL — i{Yr/Yt) sin ktL. 

This gives immediately the Bloch wavenumber 
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5. Results 


Here, we present the results of the nonlocal modeling, Bloch wave modeling and FEM simulations for 
the two-dimensional metamaterial made by Helmholtz resonators. Once the simplified nonlocal and Bloch 
wave modeling are validated by the results of the FEM simulations which are based on the solutions of 
the exact equations (8), we employ the same nonlocal modeling framework to compute the macroscopic 
acoustic properties of the three-dimensional material. For both 2D and 3D structures, the resonators are 
hlled with air as a viscothermal fluid. The fluid properties for all computations are indicated in Table 1. 
In 2D and 3D cases, the results relating to the wavenumber of the least attenuated mode and the effective 
bulk modulus of the material will be shown, versus a frequency adimensional parameter. Moreover, we 
will present a simple method allowing to obtain the 2D geometry roughly equivalent to the 3D material, 
regarding the macroscopic dynamic behavior of the material in the resonance regime of the fundamental 
mode. 

Table 1 

Fluid properties used in all computations. 


PO 

To 

CO 

U 

c 

K 

XO 

Cp 

7 

(kg/m^) 

{K) 

{m/s) 

(kg ms~^) 

(kg ms~^) 

(tPm-liC-l) 

(Pa-l) 

(j fcg-ijy-i) 


1.205 

293.5 

340.14 

1.84 X 10“® 

0.6 r) 

2.57 X 10-2 

7.17 X 10“® 

997.54 

1.4 


5.1. 2D structure filled with air 

For the geometry considered in Fig. 1 right, to perform the computations, we have set L = 1cm, S = 
0.2L, and cr = 0.015L. The functions p{u!, k) and k) are hrst determined within the approximations 

of our nonlocal modeling in section 3. Given these expressions, we know that according to nonlocal theory 
the possible wavenumbers in the medium will be the solutions of the dispersion relation (11). Solving the 
equation (11) by a Newton-Raphson scheme, we have checked that the obtained expressions for p{uj,k) 
and x~^(w, k) are such that a complex solution k{oj) to (11) exists, very close to the value ksioj) in (35). 
The frequency dependent effective density p{uj, k{uj)) = p{uj), and effective bulk modulus x“^(w, k{uj)) = 
are then obtained by putting k = k{ui) in the aforementioned excitation terms (sections 3.1 and 

3.2). 

Solving the equation (11) by the Newton-Raphson method, we varied frequency step by step, taking 
as initial value for fc(w) at a given frequency, the solution value obtained at the preceding frequency. 
Only for the starting frequency wq in the range of interest, we have chosen the value ksioJo) with a 10% 
discrepancy. 

To ascertain the validity of the modeling we have also performed direct FEM solving of the action- 
response problems, hence giving FEM evaluations of the functions p{uj,k) and x“^(w,A:). From these 
functions, the computation of the wavenumber of the least attenuated wave was performed in the same 
way as just seen, with the only difference that (for computation time reason) the initial k{uj) value at a 
given frequency was systematically taken to be ksioj) with 10% discrepancy. Finally, FEM evaluations of 
the frequency dependent effective density p{uj,k{uj)) = p{uj), and effective bulk modulus x“^(u;,fc(u;)) = 
were obtained by putting k = k{uj) in the aforementioned excitation terms. 

The FEM computations have been performed using FreeFem-l—I- [28], an open source tool solving partial 
differential equations. Adaptive meshing was employed. According to all of the calculations, the effective 
density remains practically constant and, therefore, does not play an important role in the macroscopic 
dynamics of this material. 
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Figure 3. Wavenumber (left) and bulk modulus (right) in terms of a dimensionless frequency, for the 2D structure filled 
with air. For the wavenumber, results by three calculations are compared: Bloch-wave modeling, nonlocal modeling, and 
nonlocal theory by FEM. 


We see in Fig.3, left, that the real and imaginary parts of k{uj) computed by nonlocal theory via 
Newton’s method converges exactly to the real and imaginary parts of ks which has been computed by a 
simple Bloch-wave modeling without any use of nonlocal theory. The horizontal axis is the dimensionless 
frequency koL/n, where feg = ujjcQ. The results based on the FEM simulations are also in good agreement 
with those obtained by the Bloch wave modeling and nonlocal modeling. The frequency range has been 
chosen so that it covers the resonance regime. In the same frequency range, Fig.3, right, shows the real 
and imaginary parts of the effective bulk modulus, computed by nonlocal FEM simulations and nonlocal 
modeling. Here also, we see excellent agreement between the two calculations. We notice the metamaterial 
behavior demonstrated in the real part of effective bulk modulus which becomes negative in a frequency 
range within the resonance regime. It is clear that the results by FEM computations based on the exact 
microscopic equations, can be considered more precise compared with our two modeling results in which 
we have applied simplifying approximations. As such, the good agreement between EEM results and 
others, validate the modeling framework. The discrepancies between the results based on the models and 
FEM simulations can be due in particular, to the fact that the model describes the admittance of the 
resonator W, without considering the length correction of the neck; what might generate errors in the 
calculation of the wavenumber. 

We observe here the same kind of behavior for the wavenumber and bulk modulus as it has been 
demonstrated experimentally in [23] (see Figs. 1 and 2 in that reference) for the case of the 3D material 
embedded in water. We have observed that removing the thermal effects by decreasing the coefficient of 
thermal conductivity k to a value close to zero, would have a negligible effect on the wavenumber and 
the effective bulk modulus. That is the case also for the second viscosity C) associated to losses in the 
compressional/dilatational motions in the bulk fluid. On the contrary, the material dynamics in terms 
of the macroscopic wavenumber and bulk modulus is quite sensitive to the values of the shear viscosity 
77 . In a frequency range, for instance, k^L/iT = 0.1 and 0.4, a maximum and minimum appear for the 
real part of the wavenumber. By decreasing the value of the shear viscosity, the maximum becomes 
sharper and finally diverges as the viscosity tends to zero, at the resonance frequency of the ideal fluid 
= co\/o’/[l(L — 21)(L — E — 21)], namely k^L/i: = 0.15 here; the minimum flattens and a band gap 
is created. As a matter of fact, the important feature here is the resonant behaviour of the structure 
which induces important values of the velocity in the neck, and thus also important viscous dissipation. 
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Near resonance from below, and at small enough rj, the corresponding neck flow become predominant 
and the effective wavelength is drastically reduced, leading to a so-called slow speed; but when the shear 
viscosity increases, the neck flow adjusts to a smaller value, eventually leading to the disappearance of 
the slow speed. The viscous losses also smooth out the extrema of the real and imaginary parts of the 
modulus in Fig. 3 right. Consequently, a wider frequency range of the negative real part of the bulk 
modulus is obtained by increasing the viscous losses. The thermal boundary layers near the cavity walls, 
where the fluid bulk modulus passes from adiabatic to isothermal value, mainly bring a small correction 
to the cavity spring constant (the cavity dimension is much greater than the boundary layer thickness). 
Therefore, their presence do not affect much the effective bulk modulus. 

As explained before the dynamics of the material will be very sensitive to the width of the neck, where a 
considerable part of the viscous losses takes place. Between the frequencies koL/n = 0.1 and 0.4, the ratio 
of the viscous boundary layer thickness to the width of the neck, insensibly changes from 0.35 to 0.39. We 
observed that to keep the similar behavior of the wavenumber and modulus, this ratio should remain in 
the same order, regardless of changing the scale of the material or the saturating fluid. The wavelength 
in air remains at least about 5 times larger than the periodicity L where koL/ir = 0.5, and the effective 
wavelength Ae// in the material decreases to Xeff/L ^ 8 at the resonance frequency koL/ir = 0.15, and 
to Xeff/L ~ 5 at koL/ir = 0.5. Although this structure represents a subwavelength material, and can be 
regarded in the large wavelength limit Ae// 3> L, the local theory based on the two-scale homogenization 
at order zero does not predict correctly the acoustics of this material. The origin of the failure is the 
presence of widely different length scales, allowing for resonances. 

Once the simplifying assumptions within our two modeling have been validated by the precise results 
of the FEM simulations, we can use the same modeling framework to treat the case of 3D material. 


5.2. 3D structure filled with air 

Here, the resonators are placed in a periodicity L = 1 cm, composed of a rectangular cavity of volume 
8.5 X 5 X 5 mm, a cylindrical neck I = 1 mm long and cr = 1 mm in diameter, and a main duct portion. 
The neck opens in the main square air duct with a E x E = 0.2L x 0.2L mm opening. The strategy of 
calculation to obtain the effective density, effective bulk modulus and the least attenuated wavenumber 
through nonlocal modeling and Bloch wave modeling, are the same as for 2D case in sections 3 and 4. We 
can consider that the z-axis is outward from the plane of the Fig. 2 which is regarded as a cross section 
of the 3D periodic unit. As before, Zwikker and Kosten’s plane waves are propagating and attenuating 
in the different parts of the geometry. The only change which should be applied in the 3D calculations 
with respect to 2D model, is related to the Zwikker and Kosten’s density and modulus which have been 
expressed for slits in Eq. (13). Here, we use the expressions (80) and (81) in [29] to obtain the Zwikker 
and Kosten’s density and bulk modulus for tubes of rectangular cross section (main conduit and cavity); 
for the neck (tube of circular cross section) we use the expression mentioned in [26, Appendix], and also 
in [27,29]. 

Fig. 5 left, shows the the real and imaginary parts of the frequency dependent wavenumber k{uj) 
associated with the least attenuted mode. The results based on the calculations of nonlocal modeling and 
Bloch wave modeling appear to be in perfect agreement. In Fig. 5 right, the real and imaginary parts 
of the frequency dependent bulk modulus K{uj) = x~^{^jk{uj)) are presented, according to nonlocal 
modeling. 

Between the frequencies koL/n = 0.05 and 0.3, the ratio of the viscous boundary layer to the diameter 
of the neck, changes from 0.15 to 0.06. The wavelength in air remains at least about 5 times larger than the 
periodicity L where koL/n = 0.5, and the effective wavelength in the material decreases to Xeff/L ^ 10 
at the resonance frequency koL/ir = 0.07, and to Xeff/L ^ 5 at koL/ir = 0.5. 
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2D equivalent of the 3D structure: We have performed a simple calculation to obtain the 2D 
structure made of Helmholtz resonators showing roughly the same macroscopic behavior as a 3D structure, 
in particular, in terms of the wavenumber of the least attenuated mode and the effective bulk modulus. 
We will determine the geometrical parameters of the 2D resonator, illustrated in Fig. 4, in terms of the 
parameters of a 3D resonator, in a way that it exhibits resonance at the same frequency and shows 
approximately the same dissipative character as the 3D structure does. 



Figure 4. Schematic of the 2D periodic unit of an array of Helmholtz resonator. Shown here are the geometrical parameters, 
which should be obtained in order to have a 2D equivalent of the 3D structure. 

To have roughly the same amount of both viscous and thermal losses in the main tube in 2D and 
3D, it suffices to equate the hydraulic radius (see [27]). Let be the width of the tube in 2D, and 
the side of the square tube cross section in 3D. For the hydraulic radius to be the same we take: 
= /i(^/2. In the same way, equating the hydraulic radius for the neck in 2D and 3D, gives the 
neck’s width in 2D, in terms of the diameter of the circular neck: 12. The surface of the 

cavity in 2D, 5”^ is determined in an intutive manner by assuming that the ratio of the cavity volume 
1/?-° to the tube volume in 3D is equal to the ratio of the cavity surface to the tube surface 
S'P in 2D : We will have = K F inally, the equa lity of the 

resonance frequency in 2D, uj’jf = cq \/j[K^^ and in 3D, uj^ = cq^//{ h'^V^^), results in 
the expression for the neck’s length in 2D: h'^ = 5'c^), where is the length of 

the neck in the 3D structure. 

With the dimensions of our actual 3D structure, we find, for the geometrical parameters of the 2D 
version, = 1mm, = 0.24mm, and = 8 mm. We have chosen the width = 8.5mm 
and the height = 6.25mm, so that the product of them be fixed by the value of S'^. We have 
also taken the same periodicity L for 2D, as in 3D. The complex wavenumber associated with the least 
attenuated mode, and complex effective bulk modulus of this 2D equivalent of the 3D material is depicted 
in Fig. 5. The complex wavenumber relating to the 2D and 3D geometries present an excellent agreement, 
and a very good agreement is observed regarding the real and imaginary parts of the effective bulk 
modulus of these two structures. 

We note that, if the structure with the same geometrical parameters is embedded in water, there would 
be less loss as the the viscous boundrary layer thickness is smaller compared with air. To keep the same 
dynamic behavior with water as with air, it would be necessary to very significantly decrease the width 
of the neck; at this point it should be born in mind that the complicate effect of nonlinearities would 
certainly have to be introduced. 

Note that the thermal effects in water are not important. The general thermodynamic identity 7 — 1 = 
PqTq/ pQCp, shows that the deviation of 7 = Cp/cy from unity, is a second order effect on the thermal 
expansion coefficient /3o . For a liquid, like water, jSo is very small; what implies that 7 is practically 1. In 
this case, isothermal and adiabatic bulk modulus coincide since in general K^^adiab) = iK^f^isoth)] thermal 
exchanges have virtually no effects. 
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Figure 5. Wavenumber (left) and bulk modulus (right) in terms of a dimensionless frequency, for the 3D structure filled 
with air. For the wavenumber, results by two calculations are compared: Bloch-wave modeling, nonlocal modeling. 

6. Conclusion 

Applying the Maxwellian nonlocal theory of sound propagation in porous media to a material with 
the microgeometry of the porous matrix in the form of a two or three dimensional array of Helmholtz 
resonators embedded in air, we have described precisely the metamaterial behavior of the dissipative 
medium, demonstrated by the negative real part of the effective bulk modulus in the resonance frequency 
regime. Using the homogenization method corresponding to the recently developed nonlocal theory, we 
took advantage of a plane wave modelling to obtain the effective density and bulk modulus, functions 
of both frequency and wavenumber. In this modeling we made use of Zwikker and Kosten’s equations, 
governing the pressure and velocity fields’ dynamics averaged over the cross-section of the different parts 
of the Helmholtz resonators, in order to coarse-grain them to the scale of the periodic cell containing one 
resonator. Once these two effective parameters have been determined, the corresponding least attenuated 
wavenumber of the medium could be obtained through a dispersion equation established via nonlocal 
theory. The frequency range has been chosen such that the geometrical-based resonance phenomena 
could appear. 

On the other hand, a direct analytical modelling has also been performed to obtain the least attenuated 
Bloch mode propagating in the medium, without using nonlocal theory. We have shown that the values 
of Bloch modes obtained in the direct way, match exactly those computed by the nonlocal modeling. In 
addition, the FEM numerical simulations allowing to compute the effective parameters and wavenumbers 
without any approximation, validate the results of the two modeling calculations and their simplifying 
assumptions. The nonlocal theory takes fully into account all viscous and thermal dissipation. But we have 
observed that for this material, thermal effects are negligible, while viscous effects are quite important 
to describe the material effective dynamics. We have used the same modeling framework for 3D material 
to compute the effective parameters and the wavenumber of the least attenuated mode, and performed a 
simple calculation to find a 2D equivalent of the material, showing the same macroscopic dynamics. 

Finally, the resonance induced metamaterial behavior that we have studied here, can be interpretated 
as a demonstration of the importance of considering the spatial dispersion in the medium. Higher order 
modes propagating and attenuating in this material can also be computed by the nonlocal theory and 
the subject of future research. Also, the case of the material filled with water will be analysed. 
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